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Abstract 

We consider a directed version of the Barabasi-Albert scale-free network (SFN) model with sym¬ 
metric preferential attachment with the same in- and out-degrees, and study emergence of sparsely 
synchronized rhythms for a fixed attachment degree in an inhibitory population of fast spiking 
Izhikevich interneurons. Fast sparsely synchronized rhythms with stochastic and intermittent neu¬ 
ronal discharges are found to appear for large values of J (synaptic inhibition strength) and D 
(noise intensity). For an intensive study we fix J at a sufficiently large value, and investigate the 
population states by increasing D. For small D, full synchronization with the same population- 
rhythm frequency fp and mean firing rate (MFR) /j of individual neurons occurs, while for large D 
partial synchronization with fp > {fi) {{fi): ensemble-averaged MFR) appears due to intermittent 
discharge of individual neurons; particularly, the case of fp > 4(/j) is referred to as sparse syn¬ 
chronization. For the case of partial and sparse synchronization, MFRs of individual neurons vary 
depending on their degrees. As D passes a critical value D* (which is determined by employing an 
order parameter), a transition to unsynchronization occurs due to destructive role of noise to spoil 
the pacing between sparse spikes. For D < D*, population synchronization emerges in the whole 
population because the spatial correlation length between the neuronal pairs covers the whole sys¬ 
tem. Furthermore, the degree of population synchronization is also measured in terms of two types 
of realistic statistical-mechanical measures. Only for the partial and sparse synchronization, contri¬ 
butions of individual neuronal dynamics to population synchronization change depending on their 
degrees, unlike the case of full synchronization. Consequently, dynamics of individual neurons re¬ 
veal the inhomogeneous network structure for the case of partial and sparse synchronization, which 
is in contrast to the case of statistically homogeneous random graphs and small-world networks. 
Finally, we investigate the effect of network architecture on sparse synchronization for fixed values 
of J and D in the following three cases: (1) variation in the degree of symmetric attachment (2) 
asymmetric preferential attachment of new nodes with different in- and out-degrees (3) preferential 
attachment between pre-existing nodes (without addition of new nodes). In these three cases, both 
relation between network topology (e.g., average path length and betweenness centralization) and 
sparse synchronization and contributions of individual dynamics to the sparse synchronization are 
discussed. 

PACS numbers: 87.19.1m, 87.19.1c 
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I. INTRODUCTION 


Recently, brain rhythms in health and disease have attracted much attention HIE]. Par¬ 
ticularly, we are concerned about fast sparsely synchronized brain rhythms which are related 
to diverse cognitive functions (e.g., sensory perception, feature integration, selective atten¬ 
tion, and memory formation) [3]. At the population level, synchronous small-amplitude 
fast oscillations (e.g., gamma rhythm (30-100 Hz) during awake behaving states and rapid 
eye movement sleep and sharp-wave ripple (100-200 Hz) during quiet sleep and awake im¬ 
mobility) have been observed in local held potential recordings, while at the cellular level 
individual neuronal recordings have been found to exhibit stochastic and intermittent spike 
discharges like Geiger counters HHIO]. Thus, single-cell bring activity differs distinctly from 
the population oscillatory behavior. We note that these sparsely synchronized rhythms are 
in contrast to fully synchronized rhythms where individual neurons hre regularly at the 
population frequency like the clocks. Brunei et ah developed a framework appropriate 
for description of fast sparse synchronization [IIHIS]. Under the condition of strong ex¬ 
ternal noise, suprathreshold spiking neurons discharge irregular brings as Geiger counters, 
and then the population state becomes unsynchronized. However, as the inhibitory recur¬ 
rent feedback becomes sufficiently strong, this asynchronous state may be destabilized, and 
then a synchronous population state with stochastic and intermittent individual discharges 
emerges. Thus, under the balance between strong external excitation and strong recurrent 
inhibition, fast sparse synchronization was found to occur in both random networks mHH] 
and globally-coupled networks dsnin]. 

In brain networks, architecture of synaptic connections has been found to have complex 
topology (e.g., small-worldness and scale-freeness) which is neither regular nor completely 
random [IMS]. In our recent work [26], as a complex network we employed the Watts- 
Strogatz model for small-world networks which interpolates between regular lattice with 
high clustering and random graph with short path length via rewiring [27H2^ . The Watts- 
Strogatz model may be regarded as a cluster-friendly extension of the random network by 
reconciling the six degrees of separation (small-worldness) [301 IM] with the circle of friends 
(clustering). We investigated the ebect of small-world connectivity on emergence of fast 
sparsely synchronized rhythms by varying the rewiring probability from short-range to long- 
range connection [21] • When passing a small critical value of the rewiring parameter, fast 
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sparsely synchronized population rhythms were found to emerge in small-world networks 
with predominantly local connections and rare long-range connections. We note that these 
small-world networks as well as random graphs are statistically homogeneous because their 
degree distributions show bell-shaped ones. However, brain networks have been found to 
show power-law degree distributions (i.e., scale-free property) in the rat hippocampal net¬ 
works [32H35] and the human cortical functional network [36]. Moreover, robustness against 
simulated lesions of mammalian cortical anatomical networks |371H2] has also been found 
to be most similar to that of a scale-free network (SFN) [12|. This type of SFNs are in¬ 
homogeneous ones with a few “hubs” (superconnected nodes), in contrast to statistically 
homogeneous networks such as random graphs and small-world networks |MIT5] . Many re¬ 
cent works on various subjects of neurodynamics have been done in SFNs with a few percent 
of hub neurons with an exceptionally large number of connections P61H9] . 

The main purpose of our study is to extend previous works on sparse synchronization 
in statistically homogeneous networks mHisiEni to the case of inhomogeneous SFNs with 
a few superconnected hubs. We hrst consider a directed version of the Barabasi-Albert 
SFN model with symmetric preferential attachment with the same in- and out-degrees 
= la)- [SI US] EQ] , and study emergence of sparsely synchronized rhythms 
by varying J (synaptic inhibition strength) and D (noise intensity) for a hxed attachment 
degree la in an inhibitory population of fast spiking (FS) Izhikevich interneurons [SlTIM] . 
Fast sparsely synchronized rhythms are found to appear for large values of J and D. For 
a sufficiently large hxed value of J, we make an intensive investigation of the population 
states by increasing D. For small D, full synchronization with the same population-rhythm 
frequency fp and mean bring rate (MFR) /j of individual neurons occurs. For this case, all 
the individual neurons exhibit the same behavior, independently of inhomogeneous network 
structure. As D passes a lower threshold Dth,i, a transition to partial synchronization with 
fp > ifi) {{fi)- ensemble-averaged MFR) appears due to intermittent discharge of individual 
neurons. With increasing from Dth,i, difference between fp and (fi) increases, and sparse 
synchronization with fp > 4(/j) emerges when passing a higher threshold Dth,h- For the case 
of partial and sparse synchronization, MFRs of individual neurons vary depending on their 
degrees. As D is further increased and eventually passes a critical value D*, a transition 
to unsynchronization occurs due to destructive role of noise to spoil the pacing between 
sparse spikes. The critical value D* for the transition to unsynchronization is determined 
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by employing a realistic “thermodynamic” order parameter, based on the instantaneous 
population spike rates (IPSR) It is also shown that for D < D*, population synchro¬ 
nization emerges in the whole population because the spatial correlation length between 
the neuronal pairs covers the whole system. Furthermore, the degree of the population 
synchronization is also measured in terms of two types of realistic “statistical-mechanical” 
measures, based on (1) the occupation and the pacing degrees of the spikes and (2) the 
correlations between the IPSR and the instantaneous individual spike rates |55l |56]. Only 
for the partial and sparse synchronization, contributions of individual neurons to population 
synchronization change depending on their degrees, unlike the case of full synchronization. 
Consequently, individual neuronal dynamics reveal the inhomogeneous network structure 
for the case of partial and sparse synchronization, which is in contrast to the case of sta¬ 
tistically homogeneous random graphs and small-world networks. As a next step, we also 
investigate the effect of network architecture on sparse synchronization for hxed values of 
J and D in the following three cases: (1) variation in the degree of symmetric attachment 
(2) asymmetric preferential attachment of new nodes with different in- and out-degrees (3) 
preferential attachment between pre-existing nodes (without addition of new nodes). As the 
degree la of symmetric preferential attachment in the hrst case of network architecture is 
increased, both the average path length Lp and the betweenness centralization Cb decrease, 
which results in increased efficiency of communication between nodes. Consequently, the 
degree of sparse synchronization becomes higher. On the other hand, with increasing 
the axon “wire length” of the network also increases. At an optimal degree /*, there is a 
trade-off between the population synchronization and the wiring economy, and consequently 
an optimal fast sparsely-synchronized rhythm is found to emerge at a minimal wiring cost 
in an economic SFN. As the second case of network architecture, we consider an asymmetric 
preferential attachment of new nodes with different in- and out-degrees ^ For 

this asymmetric case, we also measure Lp and Cb by varying the “asymmetry” parameter 
Ala denoting the deviation from the above symmetric case, and examine how sparse syn¬ 
chronization varies. As the magnitude |A/q,| of asymmetry parameter is increased, both Lp 
and Cb increase, which leads to decrease in efficiency of communication between nodes. As 
a result, the degree of sparse synchronization decreases. For both cases of the positive and 
the negative asymmetries with the same magnitude (e.g., AZq= 15 and -15), their values of 
Lp and Cb are nearly the same because both the inward and the outward edges are equally 
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involved in computation of Lp and Cb- However, their synchronization degrees become dif¬ 
ferent because of their distinctly different in-degree distributions affecting individual MFRs. 
In addition to the above process where preferential attachment is made to newly added 
nodes with probability a, as the third case of network architecture we also consider another 
process where preferential attachment between pre-existing nodes (without addition of new 
nodes) is made with probability (3 {a + (3 = 1). By varying /3, we also measure Lp and Cb 
and investigate the effect of this /^-process on sparse synchronization. As f3 is increased, 
communication between pre-existing neurons becomes more efficient due to decrease in both 
Lp and Cb, and hence the degree of sparse synchronization increases. For these three cases 
of network architecture, dynamics of individual neurons reveal the inhomogeneous structure 
of the SFN and hence their contributions to sparse synchronization vary depending on their 
degrees, in contrast to the case of statistically homogeneous random graphs and small-world 
networks. 

This paper is organized as follows. In Sec. we describe a directed SFN of inhibitory FS 


Izhikevich interneurons. In Sec. Ill we first investigate emergence of sparsely synchronized 
rhythms in a directed Barabasi-Albert SFN, and then the effect of network architecture 
(such as the degree of symmetric attachment, the asymmetric attachment, and the preferen¬ 
tial attachment between pre-existing nodes) on fast sparse synchronization is also studied. 


Finally, a summary is given in Section IV 


II. SCALE-FREE NETWORK OF INHIBITORY FS IZHIKEVICH INTERNEU¬ 
RONS 

We consider an SFN of N inhibitory interneurons equidistantly placed on a one¬ 
dimensional ring of radius N/2'k. Here, we employ a directed variant of the Barabasi-Albert 
SFN model, composed of two independent a— and /3—processes which are performed with 
probabilities a and (3 {a + (3 = 1), respectively [HI ESj [50]. The diagrams for these two 
processes generating an SFN are shown in Fig. 1. The a-process corresponds to a directed 
version of the Barabasi-Albert SFN model (i.e. growth and preferential directed attach¬ 
ment). For the a-process (occurring with the probability a), at each discrete time t a 
new node is added, and it has incoming (afferent) edges and outgoing (efferent) 
edges through preferential attachments with (pre-existing) source nodes and (pre- 
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existing) target nodes, as shown in Fig. 1(a). The (pre-existing) source and target nodes 
i (which are connected to the new node) are preferentially chosen depending on their out- 
degrees and in-degrees according to the attachment probabilities Usourceid-i^^^'^) 

and Utargetid'f'^^), respectively: 
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where Nt-i is the number of nodes at the time step t 
and will be referred to as symmetric and asymmetric preferential attachments, 

respectively. For the /3-process (occurring with the probability /3), there is no addition 
of new nodes (i.e., no growth), and symmetric preferential attachments with the same in- 
and out-degrees 1^)] are made between Ip pairs of (pre-existing) source and 

target nodes which are also preferentially chosen according to the attachment probabilities 
Usourceidt''^) and Utargetid^"^) of Eq. 0 , respectively, such that self-connections (i.e., loops) 
and duplicate connections (i.e., multiple edges) are excluded [see Fig. 1(b)]. Through the 
/3-process, degrees of pre-existing nodes are more intensihed. For generation of an SFN with 
N nodes, we start with the initial network at t = 0, composed of Nq = 50 nodes where the 
node 1 is connected bidirectionally to all the other nodes, but the remaining nodes (except 
the node 1) are sparsely and randomly connected with a low probability p = 0.1. Then, 
the a— and /3—processes are repeated until the total number of nodes becomes N. For our 
initial network, the node 1 will be grown as the hub with the highest degree. However, the 


results (given in Sec. Ill) are independent of the initial networks. 

As an element in our neural system, we choose the FS Izhikevich interneuron model which 
is not only biologically plausible, but also computationally efficient jSTHM] • The population 
dynamics in our SFN is governed by the following set of ordinary differential equations: 
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Here, the state of the ith neuron at a time t is characterized by two state variables: the 
membrane potential Vi and the recovery current Ui. In Eq. (|^, C is the membrane capaci¬ 
tance, Vr is the resting membrane potential, and Vt is the instantaneous threshold potential. 
After the potential reaches its apex (i.e., spike cutoff value) Vp, the membrane potential and 
the recovery variable are reset according to Eq. Q. The units of the capacitance C, the 
potential v, the current u, and the time t are pF, mV, pA, and ms, respectively. 

Unlike Hodgkin-Huxley-type conductance-based models, the Izhikevich model matches 
neuronal dynamics by tuning the parameters instead of matching neuronal electrophysiology. 
The parameters k and b are associated with the neuron’s rheobase and input resistance, a is 
the recovery time constant, c is the after-spike reset value of v, and d is the total amount of 
outward minus inward currents during the spike and affecting the after-spike behavior (i.e., 
after-spike jump value of u). Tuning these parameters, the Izhikevich neuron model may 
produce 20 of the most prominent neuro-computational features of cortical neurons [HUjM] ■ 
Here, we use the parameter values for the FS interneurons (which do not Ere postinhibitory 
rebound spikes) in the layer 5 Rat visual cortex [53|; C = 20, Vr = —55, Vt = —40, Vp = 
25, Vb = -55, k = l, a = 0.2, b = 0.025, c = -45, d = 0. 

Each Izhikevich interneuron is stimulated by using the common DC current Inc (mea¬ 
sured in units of pA) and an independent Gaussian white noise [see the 3rd and the 4th 
terms in Eq. ([^] satisfying {^i(t)) = 0 and {^i(t) ^j(t')) = 6ij 6(t — t'), where (■ ■ ■) denotes 
the ensemble average. The noise is a parametric one that randomly perturbs the strength 
of the applied current Jdc, and its intensity is controlled by using the parameter D (mea¬ 
sured in units of pA ■ ms^/^). In the absence of noise (i.e., D = 0), the Izhikevich interneuron 
exhibits a jump from a resting state to a spiking state via subcritical Hopf bifurcation for 
kDC,h = 73.7 by absorbing an unstable limit cycle born via a fold limit cycle bifurcation for 
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Idc,i = 72.8. Hence, the Izhikevich interneuron shows type-II excitability because it begins 
to hre with a non-zero frequency [371133] . As I^c is increased from lDC,hi the mean bring 
rate / increases monotonically. Throughout this paper, we consider a suprathreshold case 
oi I DC = 1500, where the membrane potential v oscillates very fast with / = 633 Hz; for 
more details, refer to Fig. 1 in [26] . 

The last term in Eq. ([^ represents the synaptic coupling of the network. Isyn,i of Eq. ([^ 
represents a synaptic current injected into the Ah neuron; gsyn,i represents the synaptic 
conductance of the Ah neuron. The synaptic connectivity is given by the connection weight 
matrix W {={wij}) where Wij = 1 if the neuron j is presynaptic to the neuron z; otherwise, 
Wij = 0. Here, the synaptic connection is modeled by using the directed SEN (explained 
in the above). Then, the in-degree of the Ah neuron, (i.e., the number of synaptic 
inputs to the neuron i) is given by '^ij■ The fraction of open synaptic ion 

channels at time t is denoted by s{t). The time course of Sj{t) of the jth neuron is given 
by a sum of delayed double-exponential functions E{t — — ti) [see Eq. (j^], where ti 

is the synaptic delay, and and Fj are the /th spike and the total number of spikes of 
the jth neuron at time t, respectively. Here, E{t) [which corresponds to contribution of a 
presynaptic spike occurring at time 0 to s{t) in the absence of synaptic delay] is controlled 
by the two synaptic time constants: synaptic rise time Tr and decay time Td, and 0(t) is the 
Heaviside step function; 0(f) = 1 for f > 0 and 0 for f < 0. For the inhibitory GABAergic 
synapse (involving the GABAa receptors), t; = 1 ms, = 0.5 ms, and = 5 ms [T3]. The 
coupling strength is controlled by the parameter J (measured in units of /rS), and Vsyn is 
the synaptic reversal potential. Here, we use Vsyn = —80 mV for the inhibitory synapse. 

Numerical integration of Eqs. ([^-(|^ is done using the Heun method [59] (with the time 
step At = 0.01 ms). For each realization of the stochastic process, we choose a random 
initial point [nj(0), Mj(0)] for the Ah (z = 1,..., A) neuron with uniform probability in the 
range of z;j(0) G (—50, —45) and zzj(O) G (10,15). 

III. EMERGENCE OF FAST SPARSELY SYNCHRONIZED RHYTHMS IN 
SCALE-FREE NETWORKS 

In this section, we study emergence of sparsely synchronized rhythms with stochastic and 
intermittent neuronal discharges by varying J (synaptic inhibition strength) and D (noise 
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intensity) in SFNs with a few superconnected hubs. Fast sparsely synchronized rhythms 
are thus found to appear for large values of J and D by employing both a thermodynamic 
order parameter and a spatial correlation function between neuronal pairs. The degree 
of population synchronization is also characterized in terms of two statistical-mechanical 
spiking and correlation measures. For this sparse synchronization, contributions of individual 
neurons to population synchronization vary depending on their degrees, and hence individual 
neuronal dynamics reveal the inhomogeneous network structure, in contrast to the case of 
statistically homogeneous random graphs and small-world networks. Furthermore, we also 
investigate the effect of network architecture on sparse synchronization for hxed J and D 
by varying (i.e., degree of symmetric preferential attachment) and AZq, (i.e., asymmetry 
parameter representing the deviation from the symmetric case) in the a-process of adding 
new nodes and the probability /3 for the /^-process of preferential attachment between (pre¬ 
existing) nodes (without addition of new nodes). 

We hrst study a directed version of the Barabasi-Albert SFN model with symmetric 
preferential attachment of = /„ = 25, composed of N inhibitory FS Izhikevich 

interneurons equidistantly placed on a one-dimensional ring of radius N/2'k |lHll5l[5n]. The 
in-degree and the out-degree of individual neurons i show power-law distributions 
with the same exponent 7 = 3.0 05], and the average number of synaptic inputs per 

neuron Ms^ (= (■ ■ ■) denotes an ensemble-average over all neurons) is 50, which is 

nearly the same as that in the small-world network of Ref. [2S]- By changing J and D, we 
investigate occurrence of population synchronized states. In computational neuroscience, an 
ensemble-averaged global potential Vg, 

(») 

^=1 

is often used for describing emergence of population synchronization. However, to directly 
obtain Vq in real experiments is very difficult. To overcome this difficulty, instead of Vg, 
we use an experimentally-obtainable IP SR (instantaneous population spike rate) which is 
often used as a collective quantity showing population behaviors [3l[IIHIS|. The IPSR is 
obtained from the raster plot of neural spikes which is a collection of spike trains of individual 
neurons. Such raster plots of spikes, where population spike synchronization may be well 
visualized, are fundamental data in experimental neuroscience. For the synchronous case, 
“stripes” (composed of spikes and representing population synchronization) are found to be 
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formed in the raster plot. Hence, for a synchronous case, an oscillating IPSR appears, while 
for an unsynchronized case the IPSR is nearly stationary. To obtain a smooth IPSR, we 
employ the kernel density estimation (kernel smoother) [50] • Each spike in the raster plot is 
convoluted (or blurred) with a kernel function Kh(t) to obtain a smooth estimate of IPSR, 
R(t): 

^ N rii 

i=l s=l 

where is the sth spiking time of the Rh neuron, rii is the total number of spikes for the 
Rh neuron, and we use a Gaussian kernel function of band width h\ 


Kh{t) = 


\/^h 




oo < t < oo. 


( 10 ) 


We hrst consider the case of H = 0. For sufficiently small J, individual interneurons hre 
too fast to be synchronized. However, as J is increased from zero, MFRs /* of individual 
interneurons decrease, and eventually when J passes a critical value J*(~ 14) a transition 
to full synchronization with the same population-rhythm frequency fp and MFR /j of in¬ 
dividual neurons occurs. Figures |^al) and |^a2) show the raster plot of spikes and the 
IPSR kernel estimate R{t) for small values of J = 100 and H = 50, respectively. Clear 
stripes are formed in the raster plot, and the corresponding IPSR kernel estimate R{t) ex¬ 
hibits large-amplitude regular oscillation with population frequency fp = 200 Hz. For this 
case, individual interneurons hre regularly with the same MFR /j which is the same as the 
population frequency fp, and hence complete full synchronization with /j = fp occurs, in¬ 
dependently of inhomogeneous network structure. However, for the sparsely synchronized 
cortical rhythms, /p : (/j) ~ 4 : 1 ((/j): ensemble-averaged MFR of individual neurons), un¬ 
like the case of full synchronization mHni. Hence, when the population frequency is much 
higher than the MFR rate of individual interneurons {fp > 4 (/*)), the synchronization will 
be referred to as sparse synchronization. For sufficiently large values of J and D, sparse syn¬ 
chronization with /p > 4 (/j) appears. Figures |^bl) and[^b2) show the raster plot of spikes 
and the IPSR kernel estimate R{t) for J = 1500 and D = 450, respectively. For this case, 
the population frequency fp of R{t) is about 147 Hz [see Fig. [^cl)], while the distribution of 
MFRs fi of individual neurons is very broad [see Fig.[^c2)] and the ensemble-averaged MFR 
ifi) (= 36 Hz) is much less than the population frequency fp. Due to this stochastic and 
intermittent discharge of individual interneurons, stripes in the raster plot become sparse 
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and smeared. Consequently, the amplitude of R{t) becomes smaller. Figure [^d) shows the 
overall state diagram in the J — D plane. As D is increased, the full synchronization for 
D = 0 evolves, depending on the values of J, and eventually desynchronization occurs when 
passing a critical value D*. Plots of fp and (/i) versus D are also shown in Figs. |^el)-[^e4) 
for J = 100, 500, 1500, and 2000. For small J [J*(— 14) < J < 173], the full synchroniza¬ 
tion for H = 0 develops directly into an unsynchronized state without any other type of 
intermediate synchronization stage because fp = /* (e.g., see the case of J = 100). However, 
for J > 173, the full synchronization for H = 0 is developed into partial synchronization 
with fp > (fi) at some lower threshold value Dthi via pitchfork-like bifurcations (e.g., see 
the cases of J = 500 1500, and 2000). With increasing J, the difference between fp and 
fi increases abruptly when passing Dth,i- For J > 1440, the partial synchronization also 
evolves into sparse synchronization with fp > 4 (fi) as D passes a higher threshold Dth,h 
(e.g., see the cases of J = 1500 and 2000), and eventually when passing a critical value D*, 
transition to unsynchronization occurs. 

For further understanding, we present explicit examples for J = 1500 which show how 
the full synchronization is evolved into an unsynchronized state as D is increased. Figures 
ial)-ia5), |^bl)-|^b5), and [^cl)-|^c5) show the raster plots, the IPSR kernel estimates 
R(t), and the inter-spike interval (ISI) histograms for D = 100, 150, 450, 600, and 800, 
respectively. For D < Dth,i{— 109), full synchronization with fp = fi occurs (e.g., see the 
case of D = 100). All the individual neurons hre regularly with the same MFR fi = 67 
Hz, which is well shown in the ISI histogram with a single peak at the global period Tq 
(~ 14.9 ms) of R{t) in Fig. |^cl). Consequently, clear stripes are formed in the raster plot 
of spikes and the IPSR kernel estimate R{t) shows large-amplitude regular oscillation with 
fp = 67 Hz [see Figs. |^al)-[^bl)]. However, when passing the lower threshold Dth,i, partial 
synchronization with fp > (fi) appears. As an example, consider the case of D = 150. In 
contrast to the case of full synchronization, the ISI histogram has multiple peaks appearing at 
multiples of the period Tq (— 9.3 ms) of R{t) [see Fig. |^c2)]. Similar skipping phenomena of 
spikings (characterized with multi-peaked ISI histograms) have also been found in networks 
of coupled inhibitory neurons in the presence of noise where noise-induced hopping from 
one cluster to another one occurs [61], in single noisy neuron models exhibiting stochastic 
resonance due to a weak periodic external force E3|, and in inhibitory networks of 
coupled subthreshold neurons showing stochastic spiking coherence [541 - I66] . “Stochastic 
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spike skipping” in coupled systems is a collective effect because it occurs due to a driving 
by a coherent ensemble-averaged synaptic current, in contrast to the single case driven 
by a weak periodic force where stochastic resonance occurs. Due to this stochastic spike 
skipping, partial occupation occurs in the stripes of the raster plot. Thus, the ensemble- 
averaged MFR (fi) (~ 46 Hz) of individual interneurons become less than the population 
frequency fp (~ 107 Hz), which results in occurrence of partial synchronization. In contrast 
to the full-synchronization case of D = 100, (fi) is decreased, while fp is increased. For 
this case of partial synchronization, the density of stripes in the raster plot becomes lower 
because smaller fraction of total neurons hre in each stripes, and the stripes become smeared, 
as shown in Fig. [^a2). Thus, both the occupation and the pacing degrees of spikes in the 
raster plot decrease, and consequently a large decrease in the amplitude of R(t) occurs 
[see Fig. [^b2)]. As D is further increased and passes the higher threshold Dth,h (— 400), 
sparse synchronization with fp > 4 (/j) appears (e.g., see the cases of D = 450 and 600). 
The interval between stripes in the raster plot becomes smaller [see Figs. [^a3)-|^a4)], and 
hence the population frequency of R{t) increases (see Figs. [^b3)-[^b4); fp= 147 and 154 
Hz for D = 450 and 600, respectively). On the other hand, the ensemble-averaged MFR 
(fi) (~36 Hz) for both cases of D = 450 and 600 is a little decreased in comparison to 
the case of D = 150, which results in decrease in density of stripes. We also note that 
multiple peaks in the ISI histogram overlap and the height of the 1st peak increases, as 
shown in Figs. |^c3)-|^c4), and hence the stripes become more and more smeared. In this 
way, both the occupation and the pacing degrees of spikes (seen in the raster plot) decrease. 
Eventually, when passing the critical value D* (~ 759), a transition to unsynchronization 
occurs. As an example of unsynchronized state, consider the case of D = 800. Multiple 
peaks in the ISI histogram become overlapped completely [see Fig. |^c5)], and hence spikes 
in the raster plot are completely scattered, as shown in Fig. |^a5). Consequently, the IPSR 
kernel estimate R{t) in Fig. [^b5) becomes nearly stationary (i.e., no population rhythm 
appears). 

In addition to the population dynamics shown in Fig. we also investigate the dynamics 
of individual neurons for J = 1500 to examine whether individual dynamics reveals the 
inhomogeneous structure of the SFN. Figures |^al)-|^a4) show the time-series of membrane 
potentials Vi of the hub neuron {i = 1 with the highest degree) and the fastest and slowest pe¬ 
ripheral neurons with low degrees [i : varying depending on D). For the full-synchronization 
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case of D = 100, all the individual neurons fire regularly with the same MFR fi (~ 67 Hz), 
as shown in Fig. |^bl), and hence complete full synchronization occurs, irrespectively of in¬ 
homogeneous structure of the SFN. However, for the partial and the sparse synchronization, 
MFRs vary depending on their degrees. For the partial synchronization of D = 150, the 
MFR /i of the hub neuron (* = 1) with highest degree is 32 Hz [which is a little less than 
the ensemble-averaged MFR (fi) (~ 46 Hz)], while the MFRs /egi and /730 of the fastest 
[i = 691) and the slowest {i = 730) peripheral neurons are 87 and 18 Hz, respectively. 
Hence, MFRs of peripheral neurons with low degrees are distributed broadly around (i.e., 
above and below) the ensemble-averaged MFR (fi) [denoted by the gray line in Fig. |^b2)], 
while MFRs of most of hub neurons with high degrees are less than (fi). As D is further 
increased and passes the higher threshold Dth,h, sparse synchronization /p > 4 (/*) appears. 
For the sparse-synchronization cases of D = 450 and 600, distributions of MFRs of individ¬ 
ual neurons become more broad when compared with that for the partial-synchronization 
case of D = 150, as shown in Figs. |^b3)-|^b4). The ensemble-averaged MFR (/j) for both 
cases of D = 450 and 600 is 36 Hz which is less than that for D = 150 because more 
fraction of neurons have lower MFRs for the case of sparse synchronization. Difference in 
the MFRs of the hub neuron {i = 1) and the fastest and the slowest peripheral neurons 
can also be easily understood in terms of the time-averaged synaptic conductance gsyn,i of 
Eq. (|^. The synaptic conductance gsyn,i of the neuron i is determined mainly by MFRs 
of pre-synaptic neurons because the fraction of open synaptic ion channels is controlled 
through the double-exponential function of spikes of pre-synaptic neurons [see Eq. Qj. If 
the MFR of a pre-synaptic neuron is fast (slow), then its contribution to gsyn,i becomes 
larger (smaller), and hence more (less) inhibition can be given to the post-synaptic neuron. 
Consequently, the MFR of the post-synaptic neuron becomes slow (fast). Figures id)- 
|^c4) show the distributions of MFRs of pre-synaptic neurons for the three cases of the hub 
neuron with i = 1 (gray region) and the fastest (solid line) and the slowest (dotted line) 
peripheral neurons. For the case of full synchronization [D = 100), all pre-synaptic neurons 
have the same MFR fi (~ 67 Hz), irrespectively of degrees of neurons. However, for the 
partial and sparse synchronization, the distribution of MFRs of pre-synaptic neurons vary 
depending on post-synaptic neurons. The fastest peripheral neuron has more fraction of 
pre-synaptic neurons with slower MFRs (as shown by the solid lines) than the hub neuron 
(gray region), and hence its time-averaged synaptic conductance becomes less than that of 
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the hub neuron. Consequently, its MFR becomes faster than that of the hub neuron. On 
the other hand, the slowest peripheral neurons have more fraction of pre-synaptic neurons 
with faster MFRs (as shown by dotted lines) than the hub neuron (gray region), and hence 
its time-averaged synaptic conductance becomes more than that of the hub neuron. As a 
result, its MFR becomes slower than that of the hub neuron. In this way, for the partial 
and sparse synchronization, individual neuronal dynamics vary depending on their degrees, 
and reveal the inhomogeneous network structure. 

As is well known, a conventional order parameter, based on the ensemble-averaged global 
potential Vg, is often used for describing transition from asynchrony to synchrony in compu¬ 
tational neuroscience [671169] . Recently, instead of Vg, we used an experimentally-obtainable 
IPSR kernel estimate i?(t), and developed a realistic order parameter, which may be ap¬ 
plicable in both the computational and the experimental neuroscience [261155]. The mean 
square deviation of R(t), 

o ^ (R(t)-my, ( 11 ) 

plays the role of an order parameter O. (Here the overbar represents the time average.) The 
order parameter may be regarded as a thermodynamic measure because it concerns just 
the macroscopic IPSR kernel estimate R(t) without any consideration between R(t) and 
microscopic individual spikes. In the thermodynamic limit of N oo, the order parameter 
O approaches a non-zero (zero) limit value for the synchronized (unsynchronized) state. 
Figure l^a) shows a plot of the order parameter versus the noise intensity D. For D < D* 
(~ 759), synchronized states exist because the order parameter O become saturated to 
a non-zero limit value for N > 3 ■ 10^. As D passes the critical value D*, a transition 
to unsynchronization occurs because the values of O tends to zero as N —)■ oo. Here we 
present two explicit examples for the synchronized and the unsynchronized states. First, 
we consider the population state for D = 700. As shown in Fig. |^bl) for N = 10^, the 
raster plot shows sparse stripes of spikes, and R(t) shows a regular oscillation, although 
there are some variations in the amplitudes. As N is increased to A^ = 10^, stripes in the 
raster plot become a little more clear, and R{t) also shows a little more regular oscillation 
[see Fig. |^b2)]. Consequently, the population state for D = 700 seems to be synchronized 
because R{t) tends to show regular oscillations as N goes to the inhnity. As a second 
example, we consider an unsynchronized case of D = 800. For N = 10^, sparse spikes are 
scattered without forming any stripes in the raster plot, and R{t) exhibits noisy fluctuations 
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with small amplitude. As N is increased to 10^, sparse spikes become more scattered, and 
consequently R{t) becomes nearly stationary, as shown in Fig. [^c2). Hence the population 
state for D = 800 seems to be unsynchronized because R{t) tends to be nearly stationary 
as N increases to the inhnity. 

We further understand the above synchronization-unsynchronization transition in terms 
of the “microscopic” dynamical cross-correlations between neuronal pairs [26] . For obtaining 
dynamical pair cross-correlations, each spike train of the ith neuron is convoluted with the 
Gaussian kernel function Kh{t) of band width h to get a smooth estimate of instantaneous 
individual spike rate (IISR) rj(f): 

rii 

= ( 12 ) 

S=1 


where is the sth spiking time of the Ah neuron, rii is the total number of spikes for the 


Ah neuron, and Kh{t) is given in Eq. (10). Then, the normalized temporal cross-correlation 
function Ci^ir) between the IISRs rj(f) and rj{t) of the (i, j) neuronal pair is given by: 


Arj(f -I- T)Arj(t) 


(13) 


where Ari{t) = ri{t) — ri{t) and the overline denotes the time average. Then, the spatial 
cross-correlation Cl (L = 1,..., N/2) between neuronal pairs separated by a spatial distance 
L is given by the average of all the temporal cross-correlations between rj(f) and 
{i = 1, ...,N) at the zero-time lag [26] : 

1 ^ 

Cl = forL = l,--- ,iV/2. (14) 

i=l 

Figure l^al) shows the plot of the spatial cross-correlation function Cl versus L for N = 10^ 
in the case of full synchronization for D = 100. The spatial cross-correlation function Cl is 
nearly non-zero constant (~ 0.97) in the whole range of L, and hence the correlation length 
Tj becomes N/2 (=500) covering the whole network (note that the maximal distance between 
neurons is N/2 because of the ring architecture on which neurons exist). Consequently, the 
whole network is composed of just one single synchronized block. For N = 10^, the flatness 
of Cl in Fig. [^bl) also extends to the whole range (L = N/2 = 5000) of the network, 
and the correlation length becomes r] = 5000, which also covers the whole network. For 
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this case oi D = 100, due to constructive role of noise favoring the pacing between sparse 
spikes, the correlation length rj seems to cover the whole network, independently of N. 
Then, the normalized correlation length fj (= ^), representing the ratio of the correlation 
length r] to the network size N (i.e., denoting the relative size of synchronized blocks when 
compared to the whole network size), has a non-zero limit value, 1/2, and consequently 
full synchronization emerges in the whole network. However, as D is further increased, the 
full synchronization breaks up due to stochastic and intermittent discharges of individual 
neurons, and then partial and sparse synchronization appears. For the cases of partial 
synchronization [D = 150) and sparse synchronization {D = 450 and 600), plots of Cl 
are shown in Figs. a2)-|^a4) for N = 10^ and in Figs. |^b2)-|^b4) for N = 10^. The 
values of Cl are also nearly non-zero constants in the whole range of L, independently of 
N. Hence, the partial and sparse synchronization appears because the correlation length r] 
covers the whole network. The degree of population synchronization may be measured in 
terms of the average spatial cross-correlation degree {Cl)l given by averaging of Cl over 
all lengths L. Figure [^c) shows the plot of {Cl)l versus D. Just after break-up of the full 
synchronization, {Cl)l drops abruptly, and then decreases slowly to zero. In contrast to the 
case of population synchronization, the spatial cross-correlation functions Cl for D = 800 
and 1000 are nearly zero for both cases of = 10^ and 10^, as shown in Figs. [^dl)-|^d2) 
and Figs. [^el)-|^e2). For theses cases, due to a destructive role of noise spoiling the pacing 
between sparse spikes, the correlation lengths r] become nearly zero, independently of iV, 
and hence no synchronization occurs in the network. 

By changing D in the whole range of population synchronization, we also measure the 
degree of population synchronization in terms of a realistic statistical-mechanical spiking 
measure Mg which was developed in our recent work [55]. As shown in Figs. |^al)-[^a4), 
population spike synchronization may be well visualized in a raster plot of spikes. For a 
synchronized case, the raster plot is composed of stripes (indicating population synchro¬ 
nization), and the density and the smearing of these stripes represent the degree of the 
population synchronization. To measure the degree of the population synchronization seen 
in the raster plot, a statistical-mechanical spiking measure M^, based on R{t), was intro¬ 
duced by considering the occupation pattern and the pacing pattern of the spikes in the 
stripes jiS] • The spiking measure Mj of the ith stripe is dehned by the product of the occu¬ 
pation degree Oi of spikes (representing the density of the ith stripe) and the pacing degree 
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Pi of spikes (denoting the smearing of the ith stripe): 


M, = O, ■ Pi. 


(15) 


The occupation degree Oi in the ith stripe is given by the fraction of spiking neurons: 


0 * = 


N 


is) 


N 


(16) 


(s) 

where N- is the number of spiking neurons in the ith stripe. For sparse synchronization, 
Oj 1, while Oi = 1 for full synchronization. The pacing degree P* of each microscopic spike 
in the zth stripe can be determined in a statistical-mechanical way by taking into account 
its contribution to the macroscopic IPSR kernel estimate R{t). Each global cycle of R{t) 
begins from its left minimum, passes the central maximum, and ends at the right minimum; 
the central maxima coincide with centers of stripes in the raster plot [see Figs. [^al)-|^a4) 
and Figs. [^bl)-[^b4)]. An instantaneous global phase <h(t) of R{t) is introduced via linear 
interpolation in the two successive subregions forming a global cycle [551 [70]; for more details, 
refer to Fig. 4 in [55]. The global phase <F(t) between the left minimum (corresponding to 
the beginning point of the Ah global cycle) and the central maximum is given by 


t -t. 


(min) 


4>(t) = 27r{t - 3/2) + TT ^ for t 


, (mm) 


<t<t 


(max) 


(j = 1,2,3,...), (17) 


and <h(t) between the central maximum and the right minimum (corresponding to the be¬ 
ginning point of the (i -t- l)th cycle) is given by 


<h(t) = 27l{i — 1) -|- TT 


t — t 


(max) 


(min) _i_(max) 


for t. 


(max) 


<t<t 


t 


i+l 


(min) 

i+l 


(* = 1,2,3, 


(18) 


where is the beginning time of the Ah global cycle (i.e., the time at which the left 

minimum of R{t) appears in the Ah global cycle) and jg the time at which the maximum 
of R{t) appears in the Ah global cycle. Then, the contribution of the fcth microscopic spike 
in the Ah stripe occurring at the time to R{t) is given by cos <Ffc, where is the global 
phase at the fcth spiking time [i.e., = 4>(t^'’^)]. A microscopic spike makes the most 

constructive (in-phase) contribution to R{t) when the corresponding global phase $a: is 27rn 
(n = 0,1, 2,...), while it makes the most destructive (anti-phase) contribution to R{t) when 
<Fj is 27r(n — 1/2). By averaging the contributions of all microscopic spikes in the Ah stripe 
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(19) 


to R(t), we obtain the pacing degree of spikes in the ith stripe: 




1 . \ 

— ^COS<hfc, 
* k=l 


where St is the total nnmber of microscopic spikes in the ith stripe. By averaging Mi of 


Eq. (15) over a snfhciently large nnmber Ng of stripes, we obtain the statistical-mechanical 
spiking measnre Mg'. 


( 20 ) 

^ 2=1 

By varying D, we follow 3 x 10^ stripes for each D and characterize popnlation synchro¬ 
nization in terms of (Oj) (average occnpation degree), (Pj) (average pacing degree), and the 
statistical-mechanical spiking measnre Mg for 11 valnes of D in the synchronized region, and 
the results are shown in Figs. [^a)-[^c). In the case of full synchronization for D < Pth,z, 
{Oi)=l and {Pi) ~ 1, which results in Mg ~ 1. However, just after break-up of the full 
synchronization, the average occupation degree (Oj) drops abruptly, because of the partial 
occupation due to stochastic spike skipping, and then it saturates to a non-zero limit value 
(~ 0.23). For the case of partial and sparse synchronization, the average pacing degree (Pj) 
also decreases monotonically to zero. Consequently, the statistical-mechanical spiking mea¬ 
sure Mg abruptly drops after break-up of the full synchronization, and then slowly decreases 
to zero, which is similar to the case of the average spatial cross-correlation degree {Cl)l 
shown in Fig. [^c). In addition to the spiking measure M^, we also characterize the pop¬ 
ulation synchronization in terms of another statistical-mechanical correlation measure M^, 
based on the cross-correlations between the IPSR R{t) and the IISRs Tiit) {i = 1,..., iV) [5B] . 
This correlation-based measure Me may also be regarded as a statistical-mechanical measure 
because it quantihes the average contribution of (microscopic) IISRs to the (macroscopic) 
IPSR. The normalized cross-correlation function C'j(r) between R{t) and rj(t) is given by 


C^{r) 


AR(t + T)Ari{t) 


( 21 ) 


where r is the time lag, AR(t) = R(t) — R(t), Ari(t) = Tiit) —ri(t), and the overline denotes 
the time average. Then, the statistical-mechanical correlation measure Me is given by the 
ensemble-average of C'j(O) at the zero-time lag |56j : 

1 ^ 

( 22 ) 

2=1 
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Figure [^d) shows the plot of Me versus D. Me — 1 for the case of full synchronization. On 
the other hand, it drops abruptly just after break-up of the full synchronization, and then 
slowly decreases to zero, which is similar to the case of Mg shown in Fig. [^c). 

For further understanding of population synchronization in Fig. we also investigate 
contributions of individual neuronal dynamics to the population synchronization. Similar 


to the population occupation, pacing, and spiking measures of Eqs. (15), (16), and (19), we 
introduce a spiking measure of the ith neuron by considering the bring and the pacing 
degrees of the spikes of the ith neuron. The bring degree F^'^\ representing the degree of 
participation of the ith neuron to the stripes in the raster plot of spikes, is given by: 


Ns 


f(*) = _L _ph) 
j=i 


(23) 


where Ng is the number stripes for averaging and denotes the participation of the ith 
neuron in the jth stripe. If the ith neuron bres in the jth stripe (i.e., the spike of the ith 
neuron participates in the jth stripe), then = 1; otherwise F^*^ = 0. The pacing degree 
of the ith neuron, denoting the degree of contributions of the spikes of the ith neuron to the 
IPSR F(t), is given by: 


5(i) 






(24) 


k=l 


where is the kth spiking time of the ith neuron {k = 1,..., is the global 

phase at and is the total number of spikes of the ith neuron. Then, the spiking 

measure of the ith neuron is given by the product of the bring and pacing degrees of 
the ith neuron: 

mW = F(*)-FW. (25) 

Figures |^al)-^cl) show plots of F^^ versus the in-degree in the case 

of the full synchronization for P = 100, respectively. The values of Fh)(= 1), p0)(~ 0.99), 
and Mi*4— 0.99) are constants, independently of the in-degrees, and hence contributions 
of individual neurons to population synchronization are the same. On the other hand, 

Ph), and vary depending on the in-degrees for the partial and sparse synchroniza¬ 
tion. The bring degrees F*-*^ of individual neurons for P = 150, 450, and 600 are shown 
in Figs. [^a2)-|^a4), respectively. Due to stochastic spike skipping of individual neurons, 
they spread around their ensemble-averaged values (F*^®)) [denoted by gray lines and cor¬ 
responding to the average occupation degree (Oj) in Fig. [^a)], as in the case of MFRs in 
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Figs.|^b2)-|^b4). Hence, of individual neurons seems to be correlated with their MFRs. 
As D is increased, the ensemble-averaged bring degree (Fh)^ decreases abruptly and then 
saturates to a lower limit value, similar to the case of {Oi) in Fig. [^a). Distributions of the 
pacing degree F^*^ and the spiking measure of individual neurons also exhibit spreads 
from their ensemble-averaged values (represented by gray lines), as shown in Figs. |^b2)- 
[^b4) and Figs. [^c2)-|^c4), respectively. With increase in D, the ensemble-averaged pacing 
degree (F*^*^), corresponding to the average pacing degree (Fj) in Fig. [^b), shows a grad¬ 
ual decrease when compared to the case of (F^)). Consequently, the ensemble-averaged 
spiking measure corresponding to the population spiking measure Mg in Fig. ^c), 

abruptly drop after break-up of the full synchronization, mainly due to sudden decrease in 
the ensemble-averaged bring degree (F^)), and then slowly decreases. With increasing D, 
the relative variances of F^), and Mi*^ from their ensemble-averaged values increase. 
For additional characterization of individual dynamics, we also introduce the correlation 
measure Mc'^ of the Ah neuron, debned by the cross-correlation C'i(O) [see Eq. (21)] between 
the IPSR R{t) and the IISR rj(f) of the Ah neuron at the zero-time lag. The “individual” 
correlation measure represents the contribution of the Ah neuron to the “population” 
correlation measure of Eq. (22). Figures ^dl)-^d4) show distributions of Mc'^ versus 
the in-degree for D = 100, 150, 450, and 600, respectively. For the case of the full 
synchronization {D = 100), Mc'^ is the same independently on the in-degrees, while for the 
cases of partial {D = 150) and sparse {D = 450 and 600) synchronization Me'* spreads 
around the ensemble-average value (denoted by gray lines). As D is increased, the 

ensemble-averaged value {Me'*) decreases, while the relative variance from {Mc^) increases, 
like the case of In this way, for the partial and sparse synchronization, contributions 

of individual dynamics to population synchronization depend on their degrees, (although 
ensemble-averages of individual measures such as pb), pb)^ and give the average occu¬ 
pation degree {Oi), pacing degree (Pi), and spiking measure Mg in the whole population,) 
and reveal the inhomogeneous structure of the SFN, in contrast to statistically homogeneous 
networks such as the random graph and the small-world network. 

From now on, we investigate the ebect of network architecture on sparse synchronization 
for bxed values of J = 1500 and D = 450 in the following three cases. As the brst case 
of network architecture, we consider the ebect of the degree P of the symmetric prefer¬ 
ential attachment (/a”^ = = /„) on sparse synchronization. Figures |^al)-^a5) and 
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Figs. [^bl)-[^b5) show the raster plots and the IPSR kernel estimates R{t) for la = 15, 20, 
25, 40, and 45, respectively. For la is less than a threshold value l^ai— population 

synchronization occurs. As an example of unsynchronization, we consider a case of la = 15 
where spikes in the raster plot are completely scattered and the IPSR kernel estimate R{t) 
becomes nearly stationary, as shown in Figs. i al) and |^bl), respectively. When pass¬ 
ing the threshold value a transition to sparse synchronization occurs. For example, for 
la = 20 stripes appear in the raster plot of spikes, and the IPSR kernel estimate R(t) shows 
regular oscillation [see Figs. |^a2) and[^b2)]. As la is further increased, the stripes in the 
raster plot become more and more dense and clear, and the IPSR kernel estimates R{t) show 
larger-amplitude regular oscillations, as shown in Figs. [^a3)-|^a5) and Figs. [^b3)-[^b5), 
respectively. Hence, as R is increased, the degree of sparse synchronization becomes better. 
For characterization of the effect of R on network topology, we also study the local property 
of the SFN in terms of the in- and out-degrees. Figures [^cl)-[^c5) show the plots of the 
out-degree versus the in-degree for R = 15, 20, 25, 40, and 45, respectively. The 
in- and out-degrees are distributed nearly symmetrically around the diagonal, and with in¬ 
creasing R they are shifted upward because of increase in the in- and out-degrees. Based on 
these degree distributions, we classify the nodes into the hub group (consisting of the head 
hub with the highest degree and the secondary hubs with higher degrees) and the peripheral 
group (composed of a majority of nodes with lower degrees). As an example, we consider 
the case of R = 25, and explain how to classify the nodes into the hub and the peripheral 
groups. For this case, the histogram for fraction of nodes versus the in-degree (which is 
also similar to that for the case of out-degree is shown in Fig. |^d). The majority of 

peripheral nodes have their degrees near the peak at = 25, while the minority of hubs 
have their degrees in the long-tail part. For convenience, we choose the threshold for 
the in-degree (denoted by the vertical dotted line in Fig.j^d) and separating the hub and the 
peripheral groups) whose fraction of nodes is 0.002 (i.e., 0.2%). Similarly, we also choose 
the threshold for the out-degree, which is the same as (Hereafter, we choose 

the thresholds and for both the in- and out-degrees whose fractions of nodes are 
0.2%). For visualization, the peripheral group is enclosed by rectangles (determined by both 
thresholds and d|^“*^) in Figs. |^cl)-^c5). The hub group (outside the rectangle) is 
composed of about 100 nodes (i.e., approximately 10% of the total neurons), where the node 
1 (denoted by the open circle) corresponds to the head hub with the highest degree and the 
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other ones will be called as secondary hubs. This kind of degree distribution is a “comet- 
shaped” one; the peripheral and the hub groups correspond to the coma (surrounding the 
nucleus) and the tail of the comet, respectively. In addition to the in- and out-degrees of 
individual nodes, we study the group properties of the SFN in terms of the average path 
length Lp and the betweenness centralization Cb by varying The average path length Lp, 
denoting typical separation between two nodes in the network, is given by the average of 
the shortest path lengths of all neuronal pairs: 

N 

where kj is the shortest path length from the node i to the node j. We note that Lp 
characterizes the global efficiency of information transfer between distant nodes. In the 
network science, centrality refers to indicators which identify the most important nodes 
within the network (i.e., the centrality indices are answers to the question “which nodes 
are most central?”). Historically hrst and conceptually simplest is the degree centrality 
(explained above), which is dehned by the number of edges of a node. This degree centrality 
represents the potentiality in communication activity. Superconnected hubs participate in 
the mainstream of information flow in the network, while peripheral nodes with a few links 
makes no active participation in the communication process. Betweenness is also another 
centrality measure of a node within the network. Betweenness centrality of the node i 
denotes the fraction of all the shortest paths between any two other nodes that pass through 
the node i [HI |72] : 

N N 

where Cjkii) is the number of shortest paths from the node j to the node k passing through 
the node i and ajk is the total number of shortest paths from the node j to the node k. 
This betweenness centrality Bi characterizes the potentiality in controlling communication 
between other nodes in the rest of the network. In our SFN, the head hub (i.e., node 1) 
with the highest degree is also found to have the maximum betweenness centrality B^ax, 
and hence the head hub has the largest load of communication traffic passing through it. 
To examine how evenly the betweenness centrality is distributed among nodes (i.e., how 
evenly the load of communication traffic is distributed among nodes), we consider the group 
betweenness centralization, representing the degree to which the maximum betweenness 
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centrality Bmax of the head hub exceeds the betweenness centrality of all the other nodes. 
The betweenness centralization is given by the sum of differences between the maximum 
betweenness centrality Bmax of the head hub and the betweenness centrality Bi of other 
node i, and normalized by dividing the sum of differences with its maximum possible value 

[7IIII2]: 


„ - (jv - i)(Ar^ - 3iv + 2) 

hb jy , max^ \Bmax Bi) . 

^SiX}_^^^^{Bmax - Bi) ^ / 


N 


(28) 


2=1 


where the maximum sum of differences in the denominator corresponds to that for the 
star network. Large Cb implies that load of communication traffic is concentrated on the 
head hub, and hence the head hub tends to become overloaded by the communication 
traffic passing through it. For this case, it becomes difficult to get efficient communication 
between nodes due to destructive interference between so many signals passing through 
the head hub ra. Figures I3<*) and [^f) show the plots of the average path length Lp 
and the betweenness centralization Cb versus respectively. With increasing both Lp 
and Cb decrease monotonically to non-zero values. Decrease in Lp implies reduction in 
intermediate mediation of nodes controlling the communication in the whole network (i.e., 
reduction in total centrality Btot given by the sum of centralities of all nodes). How the 
total betweenness Btot decreases with increase in may be seen explicitly in Fig. [^g). 
The maximum betweenness Bmax of the head hub is much more reduced than the average 
centralities of the secondary hubs and the peripheral nodes, and {B)^^^-, which leads 

to decrease in differences between Bmax of the head hub and B^ of other nodes (i.e., variation 
between centralities of nodes is reduced). Hence, as the result of increase in /„, typical 
separation between two nodes in the network becomes shorter and load of communication 
traffic becomes less concentrated on the head hub (i.e., the load is more evenly distributed 
among nodes). Consequently, as C is increased, efficiency of communication between nodes 
becomes better, which may result in the increase in the degree of sparse synchronization. 


The statistical-mechanical spiking measure Mg of Eq. (20) for the synchronization degree 
(denoted by solid circles) is shown in Fig. |^h). As C is increased, the degree of sparse 
synchronization increases and tends to become saturated. However, with increasing C, the 
network axon wiring length becomes longer due to increase in the long-range connections. 
Longer axonal connections are expensive because of material and energy costs. Hence, in 
view of dynamical efficiency we search for optimal population rhythm emerging at a minimal 
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wiring cost. We then calculate the wiring length by varying la on a ring of radius r {=N/2'k) 
where nodes are placed equidistantly. The axonal wiring length, L^w \ between the node i 
and the node j is given by the arc length between two nodes i and j on the ring: 




\j -i\ for \j-i\ < f 
N —\i — i\ for \i — > y. 


Then, the total wiring length is: 
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(29) 


(30) 


where a^- is the ij element of the adjacency matrix A of the network. The connection between 
vertices in the network is represented by its iV x adjacency matrix A (= {op}) whose 
element values are 0 or 1. If aij = 1, then an edge from the vertex i to the vertex j exists; 
otherwise no such edges exists. This adjacency matrix A corresponds to the transpose 
of the connection weight matrix W in Sec. |TT| We get a normalized wiring length by 
dividing L"”' with 1= E.ii which is the total wiring length for the 

global-coupled case: 
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Open circles in the Fig. |^h) denote the normalized wiring length Cw It increases linearly 
with respect to Iq. Hence, as la is increased, the wiring cost becomes expensive. An optimal 
rhythm may emerge through tradeoff between the synchronization degree Mg and the wiring 
cost Cw- To this end, a dynamical efficiency S is given by [ISl 1^ : 


^ Synchronization Degree (Mg) 

Normalized Wiring Length (£^) 

Figure l^i) shows plot of £ versus For la = la (= 34), an optimal rhythm is found to 
emerge at a minimal wiring cost in an economic SFN. An optimal fast sparsely synchronized 
rhythm is shown in Figs. |^jl)-|^j2). Sparse stripes appear successively in the raster plot 
of spikes. Hence, the IPSR kernel estimate R{t) shows a regular oscillation at a population 
frequency fp (~ 147 Hz), while individual neurons hre stochastically and sparsely at the 
ensemble-averaged MFR (/j) = 34 Hz. 

So far, we studied the case of symmetric attachment with = la- As the second 

case of network architecture, we consider the case of asymmetric preferential attachment 
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7^ We set 1)^’ = la + ^la and = la — ^la such that la"' + la — — 

constant, and investigate the effect of asymmetric attachment on sparse synchronization by 
varying the asymmetry parameter A/^ for = 25. For comparison, the raster plot and the 
IPSR kernel estimate R(t) for the symmetric case of Al^ = 0 (i.e., la"'^ = la^^^ = 25) are 
shown in Figs. 0a2) and0b2), respectively. Figure [lO|)al) shows the raster plot for the 
case of negative asymmetric attachment with AR = —15 [i.e., = 10 and = 40]. 

When compared with the case of the symmetric attachment, the stripes in the raster plot 
are much more smeared, while they are a little more dense. In contrast, for the case of 
positive asymmetric attachment with Al^ = 15 [i.e., la^^ = 40 and = 10], the stripes 
are less smeared but more sparse in comparison to the case of symmetric attachment, as 
shown in Fig. [Io|(a3). The amplitudes of the IPSR kernel estimates R{t) for both cases 
of Ala = —15 and 15 become smaller than that for the symmetric attachment [compare 
Figs. @bl) and [I0|[b3) with Fig. [Io|[b2)]. When the two asymmetric cases are compared, 
the amplitude of R{t) for Ala = 15 is a little larger than that for Ala = —15. In this way, 
the degree of sparse synchronization becomes reduced as the magnitude of the asymmetry 
parameter |A/q,| is increased. Depending on the sign of the asymmetry parameter A/„, the 
synchronization degree also differs, in spite of the same magnitude of AR (e.g., AR = 15 and 
-15). This difference between the cases of AR = 15 and -15 occurs due to different in-degree 
distributions affecting the synaptic inputs to individual neurons [see Eq. ([^], which will be 


l(out) 


j{in) 


l{out) 


j{in) i(out) 


explained in Fig. 11 Next, we study the effect of AR on the average path length Lp and the 
betweenness centralization Ch. Figures [T^c) and 10 d) show plots of Lp and Cb versus AR, 
respectively. Both Lp and Cb increase symmetrically with increasing |A/q|, independently of 
the sign of AR. Since both inward and outward links are involved equally in computation of 
Lp and Cb, the values of Lp and Cb for both cases of different signs but the same magnitude 
(i.e., AR and - AR) become the same, unlike the above case of population synchronization 
where only the inward synaptic inputs affect. As |A/c^| is increased, mismatching between 
the in- and out-degrees of nodes is increased, which leads to increase in Lp. This increase 
in Lp implies enhancement of intermediate mediation of nodes controlling communication 


in the network (i.e., enhancement in total betweenness Btot)- As shown in Fig. (lOKe), with 
increasing |A/a| the maximum betweenness Bmax of the head hub is much more enhanced 
than the average centralities of the secondary hubs and the peripheral nodes, and 

(^)perv which leads to increase in differences between B^ax of the head hub and Rj of 
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other nodes (i.e., variation between centralities of nodes is increased). Hence, as |A/q,| is 
increased, typical separation between two nodes in the network becomes longer and load 
of communication traffic becomes more concentrated on the head hub. Consequently, with 
increasing |A/q,|, efficiency of communication between nodes becomes worse, which may 
result in decrease in the degree of sparse synchronization. However, unlike the change in Lp 
and Cb, sparse synchronization varies depending on the sign of Ala- Figures [I0|(fl)-|l0|)f2) 
show plots of the average occupation degree (O*) and the average pacing degree (Pj) versus 
Ala- As Ala is decreased from the symmetric case (i.e., Ala = 0), (Oj) increases, while 
it decreases with increasing Ala from 0. On the other hand, with decrease in Ala from 
0, (Pi) decreases much, while it increases and tends to become saturated with increase in 
Ala from 0. As a result, the statistical-mechanical spiking measure Mg, given by taking 
into consideration both the occupation and the pacing degrees, has its peak at Ala = 0 


(i.e., symmetric case), as shown in Fig. 10 f3). Hence, Mg decreases in both positive and 
negative directions with increasing |A/q| from 0. The decreasing rate depends on the sign 
of Ala'- Mg for Ala < 0 decreases more rapidly than that for Ala > 0. For example. Mg 
for Ala = 15 is higher than that for Ala = —15. For more clear presentation, we normalize 
the occupation degree, the pacing degree, and the spiking measures by dividing them with 
their ensemble-averaged values for the symmetric case. Then, the normalized occupation 
degree (Oj), pacing degree (P*), and spiking measure Mg are shown in Fig. 10'g). As A/„ is 


decreased from 0, (Oj) increases, while (Pj) decreases much more, and hence Mg decreases. 
On the other hand, as Ala is increased from 0, (P*) increases, while (Oi) decreases much 
more, and hence Mg also decreases. Furthermore, since the variation from the symmetric 
case is larger for the case of Ala < 0, its spiking measure Mg becomes less than that for the 
positive asymmetric attachment with the same magnitude (e.g.. Mg for Ala = —15 is less 
than that for Ala = 15). 

To understand how the sparse synchronization varies differently depending on the sign 
of the asymmetry parameter Ala, we also investigate contributions of individual neuronal 
dynamics on the population synchronization. We hrst consider the effect of Ala on the degree 
distribution of nodes. Figures [ll|(al)-pT|(a3) show plots of the out-degree versus the 
in-degree for Ala = —15, 0, and 15, respectively. A majority of peripheral nodes 
with lower degrees are enclosed by rectangles, while hubs with higher degree lie outside the 
rectangles. For the case of symmetric attachment (i.e., Ala = 0), the in- and out-degrees 




are distributed nearly symmetrically around the diagonal. Hence, the in-degrees of the hubs 
and the peripheral nodes are nearly the same as the out-degrees, respectively. On the other 
hand, the degree distributions vary significantly for the case of asymmetric attachment. For 
Ala = —15, the in-degrees of peripheral nodes are less than their out-degrees, while the in¬ 
degrees of hubs are much more than their out-degrees (i.e., “popular” hubs with S> d 0 “*) 
appear). Thus, the distribution of in-degrees is broad, while the distribution of out-degrees 
is narrow (i.e., the distribution for Ala = —15 seems to be similar to that obtained through 
clockwise rotation of the symmetric distribution for Ala = 0 about a center), as shown in 
Fig. [IT|(al). In contrast, the out-degrees of peripheral nodes for Ala = 15 are less than 
their in-degrees, while the out-degrees of hubs are much more than their in-degrees (i.e., 
“social” hubs with emerge). Thus, the distribution of in-degrees is narrow, 

while the distribution of out-degrees is wide (i.e., the distribution for Ala = 15 seems to be 
similar to that obtained through counter-clockwise rotation of the symmetric distribution 
for Ala = 0 about a center), as shown in Fig. [IT|(a3). We note that individual dynamics 
vary depending on the synaptic inputs with the in-degree of Eq. (|^. Hence, the in¬ 
degree distribution affects the dynamics of individual neurons. Figures [II|(bl)-[II|(b3) show 
the power-law distributions of in-degrees for Ala = —15, 0 , and 15, respectively. As is well 
known, the exponent for Ala = 0 is 7 = 3.0 [HI HS]. On the other hand, 7 = 2.0 for 
Ala = —15 because of broad distribution, while 7 = 4.7 for Ala = 15 because of narrow 
distribution. Based on these in-degree distributions, we study MFRs of individual neurons. 
Figures [TT]( c1)-|iT|(c 3) and Figs. [TT|(dl)-|lT|(d3) show plots of MFR versus and histograms 
for fraction of neurons versus MFR for Ala = —15, 0, and 15, respectively. For the case 
of symmetric attachment (i.e., Ala = 0), the ensemble-averaged MFR (/j) [denoted by the 
horizontal gray line in Fig. [IT|[c 2 )] is approximately 36 Hz. Since the in-degree of a peripheral 
neuron is small, its pre-synaptic neurons belong to a small subset of the whole population. 
Hence, the MFRs of the peripheral neurons may change depending on the average MFR 
of pre-synaptic neurons in the small subset. If MFRs of the pre-synaptic neurons (in the 
small subset) is fast (slow) on average, then the post-synaptic peripheral neuron may receive 
more (less) synaptic inhibition, and hence its MFR becomes slow (fast). As a result, the 
MFRs of the peripheral neurons are distributed broadly around the ensemble-averaged gray 
line. The average MFR (/i)pg^j (— 38 Hz) of peripheral neurons is a little faster than the 
ensemble-averaged MFR (/*) because MFRs of the peripheral neurons are distributed a 
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little more above the horizontal gray line. On the other hand, the pre-synaptic neurons 
of a hub neuron with higher in-degree belong to a relatively larger subpopulation of the 
whole network. Since MFRs of the pre-synaptic neurons in the larger subset represent 
approximately those in the whole population, variation in the synaptic inhibitions received 
by the hub neurons is small, and hence the distribution of MFRs of the hub neurons becomes 
narrow. Moreover, since (/i)peri > (/*)> average MFR (— 25 Hz) of hub neurons 

becomes slower than the ensemble-averaged MFR (/j). Thus, MFRs of the hub neurons are 
narrowly distributed below the ensemble-averaged horizontal gray line. We then consider the 
case of the asymmetric attachment in comparison with the case of symmetric attachment. 
For Ala = —15, the in-degrees of peripheral neurons are lower, while those of hub neurons 
are much higher [compare Figs. [IT|[al) and [IT|[bl) with Figs. [IT|[a2) and [IT|[b2)]. Hence, 
the pre-synaptic neurons of a peripheral neuron belongs to a smaller subpopulation in the 
whole network. Following the same argument given in the above case of A/q, = 0, MFRs 
of the peripheral neurons are distributed around the ensemble-averaged horizontal gray line 
more broadly than those for Al^ = 0 [compare Fig. [TT|(cl) with Fig. [TT|(c2)]. As shown 
in Fig. [IT|(dl), peripheral neurons with faster MFRs appear in comparison to the case of 
Ala = 0 shown in Fig. |II|(d2), and hence the average MFR (— 50 Hz) of peripheral 

neurons becomes faster than that for Ala = 0, which also leads to increase in the ensemble- 
averaged MFR (fi) (~ 47 Hz) in the whole population, due to the majority of peripheral 
neurons. On the other hand, due to higher in-degrees, variation in the synaptic inhibitions 
received by the hub neurons becomes smaller, and hence the distribution of MFRs of hubs 
becomes more narrow. Furthermore, since (/i)perj peripheral neurons is increased, the 
average MFR {fi)hub (— ^z) of hub neurons decreases. Then, the MFRs of the hub 

neurons are more narrowly distributed much below the ensemble-averaged horizontal gray 
line [compare Fig. [IT|[cl) with Fig. [IT|[c2)]. We next consider the case of Ala = 15. For this 
case, the in-degrees of peripheral neurons are increased, while those of hub neurons are much 
decreased [compare Figs. [TT|(a3) and[TT|[b3) with Figs. [TT|[a2) an d|nfb2)] , in contrast to the 
case of Ala = —15. Hence, the pre-synaptic neurons of a peripheral neuron belongs to a little 
larger subpopulation in the whole network, and hence MFRs of the peripheral neurons are 
distributed around the ensemble-averaged horizontal gray line much narrowly than those for 
Ala = 0 [compare Fig. [IT|[c3) with Fig. [IT|[c2)]. As shown in Fig. [IT|[d3), peripheral neurons 
with slower MFRs appear in comparison to the case of Ala = 0 shown in Fig. [IT|[d2), and 
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hence the average MFR (— 29 Hz) of peripheral neurons becomes slower than that 

for Ala = 0, which also leads to decrease in the ensemble-averaged MFR (/j) (~ 28 Hz) in 
the whole population, because of the majority of peripheral neurons. Due to this narrow 
distribution of MFRs of peripheral neurons, variation in the synaptic inhibitions received 
by the hub neurons also becomes smaller, and hence the distribution of MFRs of hubs also 
becomes narrow. Moreover, since (/i)peri peripheral neurons is decreased, the average 
MFR (— 26 Hz) of hub neurons increases. Then, the MFRs of the hub neurons are 

more narrowly distributed just below the ensemble-averaged horizontal gray line [compare 
Fig. l^cS) with Fig. 002)]. 

Based on the above distributions of MFRs, we study contributions of individual dynamics 
on the sparse synchronization. Figures [ll|[el)-[ll][e3) show plots of the bring degree F*-*) of 
individual neurons versus the in-degree for Ala =-15, 0, and 15, respectively. We note 
that distributions of the bring degree of individual neurons are strongly correlated with 
their distributions of MFRs [compare Figs. [IT|(e l)-[n](e3) with Figs. [0cl)-[n](c3)]. Similar 
to the case of MFRs, F*^*^ spreads around the ensemble-averaged value (F^)) [denoted by 
gray lines and corresponding to the average occupation degree (Oi) in Fig. [Io|[fl)]. As Ala is 
increased, (F^)) decreases, which results in decrease in (Oi) in Fig. [lo|[fl). The variation of 
F*^®i about (F*^d^ ajgo decreases with increasing Ala- Distributions of the pacing degree F*^*i 
of individual neurons also show spreads from their ensemble-averaged values [represented 


by gray lines and corresponding to the average pacing degree in Fig. 10 [f2)], as shown in 
Figs.|T3fl).|I3f3), As Ala is increased, both the ensemble-averaged MFR and the variation 
decrease, and hence the ensemble-averaged pacing degree (Pbi) shows an increase, which 
also leads to increase in (Pi). Furthermore, the variation of F*-*i from (F^*i) decreases with 
increasing Ala- Figures [ll|(gl)-[ll|(g3) show plots of the individual spiking measure 
versus the in-degree for Ala =-15, 0, and 15, respectively. The value of individual 
spiking measure is determined by competition between the bring degree F*^*) and the 
pacing degree F^*^ of individual neurons, because Mi*^ is given by the product of both F*^*^ 
and F*^®^ [ see Eq. (|25[)]. For more clear comparison and presentation, we normalize the bring 
degree F^), the pacing degree and the spiking measure by dividing them with 
their ensemble-averaged values for the symmetric case of Ala = 0. The normalized bring 
degree F*^*\ pacing degree P^^\ and spiking measure are shown in Figs. [ll|[hl)-(ll|[h3). 
Figs. 3), and Figs. 01)0)3). When compared with the case of Ala = 0, 
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for Ala = “15 the normalized ensemble-averaged firing degree increases, while the 

normalized ensemble-averaged pacing degree decreases a little more. Consequently, the 
normalized ensemble-averaged spiking measure becomes less than that for Ala = 0. 

On the other hand, for Ala = 15 decreases, while (F^*^) increases only a little. As a 

result, also becomes less than that for Ala = 0. However, it is a little greater than that 

for Ala = —15 because the variation from the symmetric case of Ala = 0 is smaller for the 
case of Ala = 15. This normalized ensemble-averaged spiking measure (Afi*^) of individual 
neurons corresponds to the normalized population spiking measure Mg shown in Fig. 10 g). 


Based on the individual dynamics, it is found that the population spiking measure Mg has its 
peak value for the case of symmetric attachment due to perfect matching between the inward 
and the outward edges. As the magnitude |A/a| of the asymmetry parameter is increased 
from 0, Mg decreases in both directions because of mismatching between the inward and the 
outward edges. However, for the cases of both signs (-1-/—) with the same magnitude (e.g., 
Ala = 15 and -15) the values of Mg are different, although their network topology such as Lp 
and Cb are the same. As shown above. Mg for the case of positive asymmetric attachment 
with Ala = 15 is larger than that for the case of negative asymmetric attachment with 
Ala = —15 due to the difference in the distributions of the in-degrees. 

As the third case of network architecture, we consider the /3-process (occurring with 
the probability /3), in addition to the above a-process (which occurs with the probability a) 
{a+f3 = 1). Unlike the case of a-process, no new nodes are added, and symmetric preferential 
attachments with the same in- and out-degrees 1^)] are made between 1^ pairs 

of (pre-existing) source and target nodes which are also preferentially chosen according to the 




attachment probabilities ^source{,d!i^^"^) and of Eq. ([^, respectively, such that 

self-connections (i.e., loops) and duplicate connections (i.e., multiple edges) are excluded, 
as shown in Fig. 1(b). Here we set Iji = 5. We investigate the effect of the /3-process on 
sparse synchronization by varying /3 for the three cases of Ala =-15, 0, and 15. Figures 
il al)-[T^a3) show the raster plots of spikes for /3 =0, 0.6, and 0.8, respectively in the 
case of Ala = —15. As (3 is increased from 0, the stripes in the raster plot become more 
clear, and the IPSR kernel estimates R{t) show larger-amplitude regular oscillations, as 
shown in Figs. [I^bl)-[I^b3). Also for both cases of Ala = 0 and 15, similar effect of (3- 
process occurs in the raster plots of spikes and the IPSR kernel estimates F(t), as shown 
in Figs. [l^cl)-[l^f3). Consequently, with increasing [3 the degree of sparse synchronization 
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becomes better. For characterization of the effect of f3 on the network topology, we also 
measure the average path length Lp and the betweenness centralization by varying (3. 
Figures [T^g) and[I^h) show the plots of Lp and Ci, versus /?, respectively for the three cases 
of A/a = —15, 0, and 15. As (3 is increased, both Lp and Cb decrease monotonically for all 
three cases of A/^. As explained above, decrease in Lp leads to reduction in total centrality 
Btot (he., the sum of centralities of all nodes). How Btot decreases with increase in (3 can be 
seen explicitly in Figs. [I^il)-[I^i3) for the cases of A/^ =-15, 0, and 15, respectively. We 
note that the maximum betweenness B^ax of the head hub is much more reduced than the 
average centralities of the secondary hubs and the peripheral nodes, and for 

each case of Ala, which results in decrease in differences between B^ax of the head hub and 
Bi of other nodes (i.e., decrease in Cb)- Hence, with increasing /3, typical separation between 
two nodes in the network becomes shorter and load of communication traffic becomes less 
concentrated on the head hub. Thus, as f3 is increased, efficiency of communication between 
nodes becomes better, which may result in increase in the degree of sparse synchronization. 
Figures [l^jl)-[l^j2) show plots of the average occupation degree (Oj) and the average 
pacing degree (Pj) versus 13. As (3 is increased, at hrst (Oj) decreases for both cases of 
Ala =-15 and 0, while it increases very little for Ala = 15. Then, they seem to approach 
each other for large (3. On the other hand, (Pi) increases markedly for all the three cases 
of Ala- Consequently, the statistical-mechanical spiking measure Mg, given by taking into 
consideration both the occupation and the pacing degrees, increases monotonically mainly 


due to marked increase in (Pj) for all three cases of Ala, as shown in Fig. 12 j3). 

Finally, we investigate contributions of individual neuronal dynamics to sparse synchro¬ 
nization by varying (3 for the three cases of Ala =-15, 0, and 15. Figures |l3|;al)@a3) 
show “comet-shaped” plots of the out-degree c/(°“h versus the in-degree for (3 =0, 0.6, 
and 0.8 in the case of Ala= -15. For each 13, peripheral nodes (correspond to the coma part 
of the comet) are enclosed by the rectangle, while hubs (corresponding to the tail part of 
the comet) lie outside the rectangle and the head hub (node 1) with the highest degree is 
represented by the open circle. In the (3— process, the probability that the head hub may 
be chosen as a source and/or a target node is low because self-connections (i.e., loops) and 
duplicate connections are excluded. Hence, there is no particular change in the degree of 
the head hub, unlike the case of a—process in Fig. On the other hand, there is a marked 
increase in the degrees of some (pre-existing) peripheral nodes and secondary hubs through 
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the /3—process, which results in the immigration of some peripheral nodes into the secondary 
hub group. As a result, with increasing {3 the tail part of the comet is intensihed (i.e., the 
secondary hub group is intensihed) because the number of secondary hubs is increased. 
Then, although the number of peripheral nodes is reduced, the size of the coma part (i.e., 
the size of the rectangle enclosing the peripheral group) increases with increasing f3 because 
both the in- and the out-degrees of peripheral nodes are increased. Figures [I^dl)-[I^d3) 
also show the power-law distributions of in-degree for (3 =0, 0.6, and 0.8, respectively. 
As (3 is increased, the exponent 7 decreases, because the secondary hub group is intensihed 
(i.e. their fraction of nodes is increased) but the fraction of peripheral nodes is decreased. 


For the other two cases of A/q, = 0 and 15 (see the middle and the right panels in Fig. 13), 
as (3 is increased the distributions of in- and out-degrees evolve in a similar way, as shown in 
Figs. [I^bl)-|I^b3) and Figs. [I^c l)-[I^c3), respectively. The main ehect of /9—process is to 
intensify the secondary hub group (without particular change in the degree of the head hub) 
(i.e., the tail part of the comet-shaped distribution is intensihed with increasing f3). Hence, 
as f3 is increased the exponent for the power-law distributions of in-degree decreases [see 
Figs. [I3|[el)-[l^e3) and Figs. [I3|(fl)-|l^f3)]. These in-degree distributions for A/q, = —15, 
0, and 15 ahect the MFRs of individual neurons. Based on the change in the in-degree 
distribution in the /3—process, we study the ehect of the /3—process on the MFRs of indi¬ 
vidual neurons. Figures [T^gl)-[T^g3) show the distribution of MFRs of individual neurons 
versus the in-degree for [3 =0, 0.6, and 0.8 in the case of A/q= -15. For each [3, the 
ensemble-averaged MFR is represented by the horizontal gray line. As (3 is increased, the 
ensemble-averaged in-degree of the peripheral neurons increases, and hence the size of the 
subset of pre-synaptic neurons of a typical peripheral neuron becomes larger. Then, the 
variation of MFRs of peripheral neurons becomes reduced particularly because the part of 
higher MFRs gradually disappears. Consequently, the ensemble-averaged MFR in the whole 
population decreases a little due to the majority of peripheral neurons. With increasing (3 
the number of hubs with higher in-degrees increases, and each hub receives less inhibition 
on average because of the decreased ensemble-averaged MFR. Consequently, distribution of 
MFRs of the hubs goes upward and approaches the ensemble-averaged gray line. For the 
symmetric case of A/q = 0, distribution of MFRs of individual neurons also evolves in a 
similar way with increasing /3, as shown in Figs. [T^hl)-[T^h3). For /3 = 0 the distribution 
of MFRs of peripheral neurons becomes more narrow than that for A/q = —15 [compare 


34 



Figs. [I^hl) withjT^gl)], due to increased average in-degree of peripheral neurons (refer to 
a detailed explanation in Fig. [IT| . As j3 is increased, the average in-degree of peripheral 
neurons also increases more, and hence the variation of MFRs of peripheral neurons becomes 
more reduced. Consequently, the ensemble-averaged MFR in the whole population decreases 
very little with increasing /3, in comparison to the case of A/q, = —15. As (3 is increased, 
the distribution of MFRs of hubs goes more upward than those for A/q = —15 because the 
hubs receive less synaptic inhibition. For the case of positive asymmetric attachment with 
A/q = 15, evolution of distribution of MFRs of individual neurons also follows a similar 
way with increasing /?, as shown in Figs. [T^il)-[T^i3). For /? = 0, distribution of MFRs 
for A/q = 15 is much narrower than those for A/q = —15 and 0, and hence with increase 
in f3 the distribution of MFRs of peripheral neurons becomes reduced a little more. Then, 
the average MFR of peripheral neurons decreases a little. On the other hand, the number 
of hubs increases, their distribution of MFRs goes upward, and hence the average MFR of 
hubs increases. For this case, the effect of hubs is a little greater than that of peripheral 
neurons, and hence the ensemble-averaged MFR in the whole population becomes increased 
very little. Based on these distributions of MFRs of individual neurons, we also study 
contributions of individual dynamics to the sparse synchronization by varying (3. Figures 
[l^jl)-[l^j3) show plots of the bring degree of individual neurons versus the in-degree 
for f3 =0, 0.6, and 0.8 in the case of A/q = —15, respectively. As mentioned in Fig 
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distribution of the individual bring degree is strongly correlated with their distribution 
of MFRs. As in the case of MFRs, F^*^ also spreads around the ensemble-averaged value 
^Fh)) [denoted by the gray line and corresponding to the average occupation degree {Oi) 
(denoted by the triangles) in Fig. [I^jl)] which decreases with increasing (3. For both cases 
of increased A/q = 0 and 15, distributions of F*^*^ also evolve in a similar way with increasing 
(3, as shown in Figs. [l^kl)-p^k3) and Figs. p^ll)-[l^l3), respectively. For A/q = 0, 
decreases a little as (3 is increased, while (F^)) for A/q = 15 increases a little. These results 
of (F^*^) lead to variation of (Oi) in Fig. [l^jl). Distributions of the pacing degree F^*^ of 
individual neurons for A/q = —15, 0, and 15 also show spreads from their ensemble-averaged 
values (represented by gray lines), as shown in Figs. [I^ml)-[I^m3), Figs. [T^nl)-[T^n3), 
and Figs.[T^ol)-[T^o3), respectively. For each A/q, the variance in the distribution of MFRs 
decreases with increasing (3, and hence the ensemble-averaged pacing degree (F*^*)), corre¬ 


sponding to the average pacing degree (Fj) in Fig. 12 j2), increases monotonically as (3 is 
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increased, while the variation of from tends to decrease. Figures [l^pl)-[l^p3), 

Figs. [I^ql)-[I^q3), and Figs. [I^rl)-p!^r3) show plots of the individual spiking 


measure 


Ms versus the in-degree for Ala =-15, 0, and 15, respectively. For all three cases of 
Ala, with increasing (5 (P^d^ increases markedly, in comparison to (pd)). Consequently, for 
each Ala the ensemble-averaged spiking measure corresponding to the population 


spiking measure Mg in Fig. 12 j3), increases as [3 is increased. 

In the way explained above, for all the three cases of effect of network architecture on 
sparse synchronization, contributions of individual neuronal dynamics on population syn¬ 
chronization vary depending on the in-degrees, although their ensemble-averaged individual 
measures, (P(d), (pW), and (mF), give the population measures such as (Oj), {Pi), and 
Mg. Consequently, dynamics of individual neurons for the case of sparse synchronization 
reveal the inhomogeneous structure of the SFN, in contrast to statistically homogeneous 
random and small-world networks. 


IV. SUMMARY 

In order to extend previous works on sparse synchronization in statistically homogeneous 
networks such as random graphs and small-world networks [TTf[T6| l26] to the case of inhomo¬ 
geneous networks, we have investigated emergence of sparsely synchronized brain rhythms in 
the directed version of the Barabasi-Albert SFN model with symmetric preferential attach¬ 
ment with the same in- and out-degrees {la^^ = la^^^ = la = 25). Fast sparsely synchronized 
rhythms with stochastic and intermittent neuronal discharges have been found to emerge for 
large values of J and D. We have made an intensive investigation of population states by 
varying D for a fixed value of J = 1500. For small D, fully synchronized rhythms with the 
same fp (population-rhythm frequency) and fi (MFR of individual neurons) appear. For this 
case of full synchronization, all the individual neurons exhibit the same oscillatory behaviors, 
independently of inhomogeneous network structure. However, as D passes a lower threshold 
Dth,i{— 109), partial synchronization with fp > (fi) (ensemble-averaged MFR of individ¬ 
ual neurons) occurs due to intermittent discharge of individual neurons. As D is increased 
from Dth,i, difference between fp and (/*) increases, and when passing a higher threshold 
Dth,h{— 400), sparsely synchronized rhythms with fp > 4(/j) appear. Unlike the case of 
full synchronization, MFRs of individual neurons vary depending on their in-degrees for the 
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case of partial and sparse synchronization. As D is further increased and eventually passes a 
critical value Z)*(~ 759), a transition to unsynchronization occurs due to destructive role of 
noise to spoil the pacing between sparse spikes. The critical value D* has been determined 
through calculation of the thermodynamic order parameter O. For D < D*, population 
synchronization has been found to emerge because the spatial correlation length between 
the neuronal pairs covers the whole system. Moreover, the degree of population synchro¬ 
nization has also been measured in terms of two types of statistical-mechanical spiking and 
correlation measures. Unlike the case of full synchronization, individual neuronal dynamics 
vary depending on their in-degrees and reveal the inhomogeneous network structure for the 
case of partial and sparse synchronization, in contrast to the case of statistically homoge¬ 
neous random graphs and small-world networks. As a next step, we have also investigated 
the effect of network architecture on sparse synchronization for hxed values of J = 1500 and 
D = 450 in the following three cases: (1) variation in the degree of symmetric attachment 
(2) asymmetric preferential attachment of new nodes with different in- and out-degrees (3) 
preferential attachment between pre-existing nodes (without addition of new nodes). As the 
degree la of symmetric preferential attachment is increased, both the average path length 
Lp and the betweenness centralization Cb have been found to decrease. Hence, typical sepa¬ 
ration between two nodes in the network becomes shorter and load of communication traffic 
becomes less concentrated on the head hub. Consequently, with increasing la the degree 
of sparse synchronization has been found to become higher due to increased efficiency of 
communication between nodes. On the other hand, the normalized axon wire length of 
the network also increases. Through a trade-off between the population synchronization and 
the wiring economy, an optimal sparsely-synchronized rhythm has been found to emerge at 
a minimal wiring cost in an economic SFN with an optimal degree la{— 34). As the second 
case of network architecture, we have also considered the case of asymmetric preferential 
attachment of new nodes with different in- and out-degrees (i.e., ^ la^^^)- For this case, 

we have also measured Lp and Cb by varying the asymmetry parameter Ala denoting the 
deviation from the above symmetric case, and investigated how the sparse synchronization 
changes. As the magnitude \Ala\ of the asymmetry parameter increases, both Lp and Cb 
have been found to increase symmetrically, independently of the sign of Ala- Hence, with 
increasing |A/q,| typical separation between two nodes in the network becomes longer and 
load of communication traffic becomes more concentrated on the head hub, due to increased 
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mismatching between the inward and outward edges. Consequently, as \Ala\ is increased 
the degree of sparse synchronization has been found to become lower due to decreased ef- 
hciency of communication between nodes. For both cases of the positive and the negative 
asymmetries with the same magnitude (e.g., Al^ = 15 and -15), the values of Lp and Cb 
have been found to be nearly the same, because both inward and outward edges are involved 
equally in computation of Lp and Cb- However, their degrees of sparse synchronization have 
been found to become different due to their distinctly different in-degree distributions af¬ 
fecting individual MFRs. In addition to the above a-process where preferential attachment 
is made to newly added nodes with probability a, we have also considered another /3-process 
where preferential attachment between pre-existing nodes (without addition of new nodes) 
is made with probability /3 {a + /3 = 1). By varying the probability /3, we have also mea¬ 
sured Lp and Cb and investigated the effect of this /3-process on sparse synchronization. As 
/3 is increased, communication between pre-existing neurons becomes more efficient due to 
decrease in both Lp and Cb, and consequently the degree of sparse synchronization has been 
found to increase. For these three cases of network architecture, contributions of individual 
neuronal dynamics on the sparse synchronization were also characterized in terms of their 
MFRs, bring degrees, pacing degrees, and spiking measures. It has thus been found that the 
dynamics of individual neurons vary depending on their in-degrees and reveal the inhomo¬ 
geneous structure of the SFN, in contrast to the statistically homogeneous networks such as 
the random graph and the small-world network. Finally, we expect that our results might 
provide important insights on emergence of fast sparsely synchronized rhythms, associated 
with diverse cognitive functions such as sensory perception, feature integration, selective 
attention, and memory formation, in real brain networks with scale-free property. 
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(a) a-process (b) p-process 



FIG. 1: Diagrams of two processes generating a directed SFN. (a) Diagram of the a-process of 
adding a new node (denoted by a gray circle) with preferential attachment of inward and 
outward edges, (b) Diagram of the /?—process of the preferential attachment between Ijs pairs of 
(pre-existing) source and target nodes without adding a new node. The open circles with labels 
“S'” and “T” represent the (pre-existing) source and target nodes, respectively. 
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FIG. 2: State diagram in the J — D plane for I^c = 1500 in the directed SFN of a = 1 (i.e., 
(3 = 0) and = /a = 25 (i.e., symmetric preferential attachment), (al) Raster plot of 

spikes and (a2) plot of the IPSR kernel estimate R{t) for the full synchronization when J = 100 
and D = 50. (bl) Raster plot of spikes and (b2) plot of the IPSR kernel estimate R{t) for the 
sparse synchronization when J = 1500 and D = 450. The band width of the Gaussian kernel 
estimate for the IPSR R{t) is 1 ms. (cl) One-sided power spectrum of AR{t) [= R{t) — R{t)] 
(the overbar represents the time average) with mean-squared amplitude normalization and (c2) 
distribution of mean hring rates (MFRs) of individual neurons for J = 1500 and D = 450. Power 
spectrum is obtained from 2^® (=65536) data points. Averaging time for the MFR is 10^ ms 
and the bin size for the histogram is 3 Hz. (d) State diagram in the J — D plane. For the full 
synchronization, the individual MFR /j is the same as the population frequency fp, while for the 
partial and sparse synchronization, the ensemble-averaged MFR (/j) is less than fp. Particularly, 
the case of fp > 4(/j) is referred to as the sparse synchronization. Plots of fp and (/j) versus for 
J = (el) 100, (e2) 500, (e3) 1500, and (e4) 2000. Here, the circles and the crosses denote fp and 
(/j), respectively. 
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FIG. 3: Fast sparse synchronization for = 1500 and J = 1500 in the directed SFN of a = 1 
(i.e., /3 = 0) and = la = 25 (i.e., symmetric preferential attachment). Raster plots of 

spikes in (al)-(a5), plots of the IPSR kernel estimate R{t) in (bl)-(b5), and the ISI histograms in 
(cl)-(c5) for various values of = 100, 150, 450, 600, and 800; vertical dotted lines denote integer 
multiples of the global period Tq of R{t) 9.3 ms in (c2), 6.8 ms in (c3), and 6.5 ms in (c4)]. 
The band width of the Gaussian kernel estimate for the IPSR R{t) is 1 ms. Each ISI histogram is 
composed of 5 x 10^ ISIs and the bin size for the histogram is 0.5 ms. 
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FIG. 4: Individual neuronal dynamics for I^c = 1500 and J = 1500 in the directed SFN of a = 1 
(i.e., /? = 0) and = /^ = 25 (i.e., symmetric preferential attachment). (al)-(a4) Time 

series of the membrane potentials of the hub neuron (i = 1) with highest degree and the fastest and 
the slowest peripheral neurons with low degrees for various values of H = 100, 150, 450, and 600. 
(bl)-(b4) Plots of MFRs of individual neurons versus the in-degree for D = 100, 150, 450, and 
600. Averaging time for the MFR is 10^ ms. (cl)-(c4) Histograms for the MFRs of pre-synaptic 
neurons for the hub neuron and the fastest and the slowest peripheral neurons with low degrees for 
D = 100, 150, 450, and 600. Horizontal gray lines in (b2)-(b4) denote ensemble-averaged MFRs [~ 
46 Hz in (b2) and 36 Hz in (b3)-(b4)]. Gray regions, solid lines, and dotted lines in (c2)-(c4) denote 
the histograms for the MFRs of the pre-synaptic neurons for the hub neuron and the fastest and 
the slowest peripheral neurons with low degrees, respectively. Histograms in (cl)-(c4) are obtained 
from 30 realizations and the bin size for the histogram is 3 Hz. 


45 












(b2) D=700, N=10" 


(b1) D=700, N=10' 



N=10'’ 1000 1025 1050 1000 1025 1050 

t (ms) t (ms) 


FIG. 5: Synchronization-unsynchronization transition for Ijyc = 1500 and J = 1500 in the directed 
SFN of a = 1 (i.e., /3 = 0) and = la = 25 (i.e., symmetric preferential attachment), 

(a) Plots of the thermodynamic order parameter versus D. Averaging time for the thermodynamic 
order parameter is 3 x 10^ ms. Synchronized state for D = 700; raster plots of spikes and plots of 
the IPSR kernel estimate R{t) for N = (bl) 10^ and (b2) 10^. Unsynchronized state for D = 800: 
raster plots of spikes and plots of the IPSR kernel estimate R{t) for N = (cl) 10^ and (c2) 10^. 
The band width of the Gaussian kernel estimate for the IPSR is 1 ms. 
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FIG. 6: Characterization of synchronization-unsynchronization transition in terms of spatial cross¬ 
correlations for Idc = 1500 and J = 1500 in the directed SFN of a = 1 (i.e., (3 = 0) and 
[On) _ = 25 (i.e., symmetric preferential attachment). Plots of the spatial correlation 

function Cl between neuronal pairs versus spatial distance L for the synchronized cases of various 
values of D =100, 150, 450, and 600 when N = (al)-(a4) 10^ and (bl)-(b4) 10^. (c) Plot of the 
average spatial-correlation degree {Cl)l versus D. Plots of the spatial correlation function Cl 
versus L for the unsynchronized cases oi D = 800 and 1000 when (dl)-(d2) N = 10^ and (el)-(e2) 
N = 10^. The number of data used for the calculation of each temporal cross-correlation function 
Cijir) (the values at zero-time lag (r = 0) are used for calculation of Cl) is 2 x 10^. 
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FIG. 7: Characterization of population synchronization for /£)C' = 1500 and J = 1500 in the 
directed SFN of a = 1 (i.e., /3 = 0) and = G = 25 (i.e., symmetric preferential 

attachment). Plots of (a) the average occupation degree (O*), (b) the average pacing degree (Pj), 
and (c) the statistical-mechanical spiking measure Mg versus D. (Oj), {Pi), and Mg are obtained 
by following 3 x 10^ stripes in the raster plot of spikes, (d) Plot of the statistical-mechanical 
correlation measure Me, based on temporal cross-correlations between the IPSR R{t) and IISRs 
ri{t) of individual neurons versus D. The number of data used for the calculation of temporal 
cross-correlation function for each P is 2 x 10^. 
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FIG. 8: Contributions of individual dynamics to the population synchronization for I^c = 1500 
and J = 1500 in the directed SFN of a = 1 (i.e., /3 = 0) and = la = 25 (i.e., symmetric 

preferential attachment). Plots of (a) the firing degree (b) the pacing degree and (c) the 
spiking measure of individual neurons versus the in-degree for various values of D =100, 
150, 450, and 600. FW, and are obtained by following 3 x 10^ stripes in the raster 
plot of spikes, (d) Plots of cross-correlations between IPSR R{t) and IISRs ri{t) of individual 
neurons versus the in-degree for D =100, 150, 450, and 600. The number of data used for 
the calculation of temporal cross-correlation function for each is 2 x 10^. Horizontal gray lines 
represent ensemble-averaged values. 
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FIG. 9: Effect of symmetric attachment degree la on the sparse synchronization and economic 
SFN for Idc = 1500, J = 1500, and D = 450 in the directed SFN of a = 1 (i.e., /3 = 0). Raster 
plots of spikes in (al)-(a5) and plots of the IPSR kernel estimate R{t) in (bl)-(b5) for various values 
of symmetric attachment degree la- The band width of the Gaussian kernel estimate for the IPSR 
R{t) is 1 ms. Plots of the out-degree versus the in-degree for la = (cl) 15, (c2) 20, (c3) 
25, (c4) 40, and (c5) 45. Peripheral groups are enclosed by rectangles, while hubs lie outside the 
rectangles. The head hub with the highest degree is represented by the open circle, (d) Histogram 
for fraction of nodes versus the in-degree for la = 25. This histogram is obtained through 
30 realizations and the bin size for the histogram is 1. The vertical line represents a threshold 
for whose fraction of nodes is 0.002 (i.e., 0.2%). Plots of (e) average path length Lp and (f) 
betweenness centralization Cb versus la- (g) Plots of the maximum betweenness centrality Bmax-, the 
average betweenness centrality of secondary hubs, and the average betweenness centrality 

{B)peri of peripheral nodes versus la- Here, (• • •)^ represents an average over 30 realizations, (h) 
Plots of statistical-mechanical spiking measure Mg and normalized wiring length Cw versus la- (i) 
Dynamical efficiency S versus la- The values of Mg, and £ at an optimal value of la = 34 are 
denoted by the symbol Optimally fast synchronized rhythm for la = la- (jl) raster plot of 
neural spikes and (j2) plot of the IPSR kernel estimate R{t). 
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FIG. 10; Effect of asymmetric attachment on the sparse synchronization for l£>c = 1500, J = 1500, 
and D = 450 in the directed SFN of a = 1 (i.e., /? = 0), = 25 + Al^, and = 25 — AIq. 

{Ala- asymmetry parameter). Raster plots of spikes in (al)-(a3) and plots of the IPSR kernel 
estimate R{t) in (bl)-(b3) for various values of Ala= -15, 0, and 15. The band width of the 
Gaussian kernel estimate for the IPSR R{t) is 1 ms. Plots of (c) average path length Lp and 
(d) betweenness centralization Cb versus Ala- (e) Plots of the maximum betweenness centrality 
Bmax, the average betweenness centrality of secondary hubs, and the average betweenness 

centrality of peripheral nodes. Here, (• • • )^ represents an average over 30 realizations. 

Plots of (fl) the average occupation degree (O*), (f2) the average pacing degree (Pi), and (f3) the 
statistical-mechanical spiking measure Mg versus Ala- {Oi), {Pi), and Mg are obtained by following 
3 X 10^ stripes in the raster plot of spikes, (g) Plots of the normalized average occupation degree 
{Oi), the normalized average pacing degree {Pi), and the normalized statistical-mechanical spiking 
measure Mg- Normalizations of {Oi), {Pi), and Mg are done by dividing them with the values for 
the case of Ala = 0. 51 











FIG. 11: Distinct differences in individual neuronal dynamics for the case of asymmetric attachment 
when Idc = 1500, J = 1500, and D = 450 in the directed SFN of a = 1 (i.e., j3 = 0), = 

25+AZq,, and = 25 — AZ^ (AZq, = -15, 0, and 15). (al)-(a3) Plots of the out-degree versus 
the in-degree cZ^*”) and (bl)-(b3) power-law in-degree distributions with different exponents for the 
cases of AZq, = -15, 0, and 15. The fractions of nodes are 0.2% at the thresholds and 
of the in- and out-degrees which determine the rectangle enclosing the peripheral group. (cl)-(c3) 
Plots of MFRs versus the in-degree and (dl)-(d3) histograms for fraction of neurons versus 
MFR for the cases of AZq = -15, 0, and 15. Plots of (el)-(e3) firing degree (fl)-(f3) pacing 
degree and (gl)-(g3) spiking measures for individual neurons versus the in-degree 
and plots of (hl)-(h3) normalized firing degree (il)-(i3) normalized pacing degree and 
(jl)-(j3) normalized spiking measure for the cases of AZq = -15, 0, and 15. Normalizations 
of f(*), and are done by dividing them with the values for the case of AZq = 0. Gray 
horizontal lines represent the ensemble-averaged values. 
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FIG. 12: Effect of /3-process on sparse synchronization for I^c = 1500, J = 1500, and D = 450 in 
the directed SFN of /3 = 0, 0.6, and 0.8 (i.e., a = 1 — /3). For the a—process, = 25 + Ala and 
_ 25 _ [Ala = -15, 0, and 15), while for the /3—process = 5. Effect of 

/3-process for Ala = —15: (al)-(a3) Raster plots of spikes and (bl)-(b3) plots of the IPSR kernel 
estimate R{t) for the cases of /3 = 0, 0.6, and 0.8. Effect of /3-process for Ala = 0: (cl)-(c3) Raster 
plots of spikes and (dl)-(d3) plots of the IPSR kernel estimate R{t) for the cases of /3 = 0, 0.6, 
and 0.8. Effect of /3-process for Ala = 15: (el)-(e3) Raster plots of spikes and (fl)-(f3) plots of the 
IPSR kernel estimate R{t) for the cases of /3 = 0, 0.6, and 0.8. Plots of (g) average path length Lp 
and (h) betweenness centralization Ch versus /3 and (il)-(i3) Plots of the maximum betweenness 
centrality Bmax, the average betweenness centrality of secondary hubs, and the average 

betweenness centrality of peripheral nodes versus 13 for the cases of Ala =-15, 0. and 15. 

Here, (•••),, represents an average over 30 realizations. Plots of (jl) average occupation degree 
(Oj), (j2) average pacing degree (Pi), and (j3) statistical-mechanical spiking measure Mg versus /3 
for the three cases of Ala = -15, 0, and 15. 53 
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FIG. 13: Distinctly different individual dynamics in the /3-process for I^c = 1500, J = 1500, 
and D = 450 in the directed SFN of /3 = 0, 0.6, and 0.8 (i.e., a = 1 — /3). For the a—process, 
= 25 + Ala and = 25 — Ala (Ala = -15, 0, and 15), while for the /3—process = 
l{^t) = Plots of the out-degree versus the in-degree for Ala= (al)-(a3) -15, 

(bl)-(b3) 0, and (cl)-(c3) 15. The fractions of nodes are 0.2% at the thresholds and of 
the in- and out-degrees which determine the rectangle enclosing the peripheral group. Power-law 
in-degree distributions for Ala= (dl)-(d3) -15, (el)-(e3) 0, and (fl)-(f3) 15. Plots of MFRs versus 
the in-degree for Ala= (gl)-(g3) -15, (hl)-(h3) 0, and (il)-(i3) 15. Plots of firing degree 
for individual neurons versus the in-degree d*^®®®^ for Ala= (jl)-(j3) -15, (kl)-(k3) 0, and (11)-(13) 15. 
Plots of pacing degree for individual neurons versus the in-degree d^®®®) for Ala= (ml)-(m3) 
-15, (nl)-(n3) 0, and (ol)-(o3) 15. Plots of spiking measure for individual neurons versus the 
in-degree d*^®®®) for Ala= (pl)-(p3) -15, (ql)-(q3) 0, and (rl)-(r3) 15. 
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